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(57) Abstract: The present invention is an incremental umbrella sampling method to 
improve the performance of established sampling methods. It is sampling the state 
space by iteratively generating states xU , and their weighting factors represented by 
Formula (a) by fitting the sampling distribution function p/x) of the next iteration to 
at least one weighted property of the already sampled states. This means that p/x) is 
fitted to the product represented by Formula (b), in which Formula (a) is the weight- 
ing factor and 0( Xti ) is a function respectively a property of the states x itt . The number 
of states Xij and the number of weighting factors (see Formula (a)) is incremented 
with each iteration. In order to have a consistent set of weighting factors (see For- 
mula (a)), the weighting factors are recalculated in each iteration for all, respectively 
for a set of selected, states. By fitting p/x) in the state space it is possible to use 
all the information of Formula (a) and 0(x itt ) for the states x h , generated so far. The 
fitting step allows to use different fitting strategies. For example the fitting can bias 
the sampling away from areas where intensive sampling has been done in the pre- 
ceding iterations, or the sampling can be directed along local gradients respectively 
towards local minima or maxima of one or several weighted properties. In each of 
the iterations, the sampling distribution function is fitted in a way to improve the 
overall sampling of the state space. The method supports multi -objective optimisa- 
tions. State space integrals can be solved. It reduces the probability that the system 
is trapped. The invention is general. It can be used with different sampling methods, 
in particular with Monte Carlo sampling, Metropolis Monte Carlo sampling, or dy- 
namic simulations. It can be combined with the concepts of simulated annealing and 
multicanonical sampling. It provides a general framework that can be adapted to the 
system and the observables of interest. 
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Sampling method 

Field of the Invention 

Identifying one or a few states with desired properties, and averaging properties over all states 
of a system, can provide solutions to many important problems including determinations of 
binding free energies of drug candidates, protein structure predictions, analysis of experimental 
data, determinations of the value at risk of financial portfolios, or option pricing. Sampling 
algorithms such as Monte Carlo sampling, Metropolis Monte Carlo sampling, or dynamic 
simulations are general methods to address complex optimisation problems and to solve state 
space integrals. The problem with sampling algorithms is that the number of states required to 
obtain converged results is often large, and sometimes prohibitively large. The systems studied 
with sampling algorithms have, in general, a large number of states x with probabilities or 
probability densities that are known at least up to a normalisation constant. Examples are 
molecular systems, where x may denote the co-ordinates of the atoms, or financial systems, 
where x may be a vector of risk factors such as interest rates or stock indices. 

For such systems, the average (Ok) of a property Ok(x) can be calculated with Eq. (1), if there 
is a continuum of states, /.e., 

(a)=jXx)a(jc)ix (D 

In Eq. (1), the integral is over the entire state space, p(x) is the probability density that the 

system is in state x, and the index k may be used to distinguish different properties. Similar, if 

the states are discrete, sums 

(0')=ZpMxj) (2) 
J 

over all states must be calculated where pj is the probability of the system to be in state xj 
Generalisations to other cases are possible. 

An example is the estimation of the average distance between two atoms of a molecular 
system. The average distance (r) is equal to the integral 

(R)=$dxp(xXx) (3) 
over all conformations of the system. In Eq. (3), x is the vector of all co-ordinates of the system, 
f{x) denotes the probability distribution of the system, i.e., p(x) is< equal to the probability 
density that the system has co-ordinates x, and R(x) is the distance between the two atoms of 
interest as a function of the co-ordinates of the system. 
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The identification of some states with desired properties is a special case of a state space 
integral. E.g., the identification of the state jc™ that corresponds to the unique global minimum 
of some function f(x) that has a unique global minimum is equal to solving the integral 

Xrrm = tyxs{f(^mm{f(xl^ . (4) 
5 In Eq. (4), min(/(^/)) denotes the unique global minimum of the function f(x). 

Prior Art 

For many problems of interest, the integrals (Eq. 1) or sums (Eq. 2) cannot be solved 
0 analytically, and sampling algorithms have to be used to estimate the properties. Sampling 
algorithms produce a sample of states such that estimates of the integral of Eq. (1) or the sum of 
Eq. (2) are obtained as a weighted sum over the states of the sample, i.e., 

&=£>^(*<)- (5) 

In Eq. (5), Ok denotes the estimate of the average (Ok) of the property Ok(x), pij is the weight 

5 of the state xu , the sum is over the N§ states produced by the sampling algorithm, the index t 
distinguishes different states, the index k can be used to distinguish different properties, and the 
index / can be used to distinguish different runs of the sampling algorithm. 

Established sampling methods include systematic search methods and random sampling 
o methods. Systematic methods, evaluate the properties of the system for a set of states 
determined in advance by a deterministic algorithm. Each state produced represents a subset of 
all states of the system. The weight used in Eq. (5) for state xr.r is proportional to the size of 
the subset represented by the state xv . Systematic search is often used to estimate low 
dimensional integrals, and improvements of the simple method just described are known and 
5 applied routinely, e.g., trapezoidal interpolation. However, there are cases where prohibitively 
large sets of states are required to obtain converged estimates from systematic searches, e.g., 
for systems with a large number of degrees of freedom. 

Random sampling methods generate random states xu distributed according to a sampling 
) distribution function. For each state considered in run /, the weight put for Eq. (5) can be 
determined, e.g., by Maximum Likelihood considerations (Christian Bartels "Analysing biased 
Monte Carlo and molecular dynamics simulations", Chem. Phys. Lett. 331 (2000) 446-454; see 
also M. Souaille and B. Roux, "Extension to the Weighted Histogram Analysis Method: 
Combining Umbrella Sampling with Free Energy Calculations", CompuL Phys. Comm. 135, 
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(2001) 40-57; Charles J. Geyer, "Estimating normalizing constants and reweighting mixtures in 
Markov chain Monte Carlo", Tech. Rep. 568r (1994), School of Statistics, University of 
Minnesota; Art B. Owen and Yi Zhou, "Safe and effective importance sampling", Journal of the 
American Statistical Association 95 (2000) 135-; A.M. Ferrenberg and R.H. Swendsen, 
"Optimized Monte Carlo Data Analysis", Phys. Rev. Lett 63, (1989) 1 195-). The weights depend 
on the sampling distribution functions used. Different random sampling methods such as Monte 
Carlo, Metropolis Monte Carlo, or Langevin Molecular Dynamics differ in the way the states are 
generated. In Monte Carlo sampling - for example according to US 6 061 662 - states are 
generated independently of each other. Metropolis Monte Carlo generates new states by 
perturbing existing states and accepting or rejecting them such that the states are sampled 
according to a given sampling distribution function. 

With the mentioned sampling algorithms, it is in principle possible to estimate properties of 
arbitrary systems. However, in many cases, the size of the sample required to obtain converged 
estimates is large, and calculation times get prohibitively long. The convergence of the sampling 
algorithm depends critically'on the sampling distribution function, which is used. 

For many systems and properties of interest, the best choice of the sampling distribution 
function is the distribution function of the system p(x). This choice of the sampling distribution 
function is often optimal to estimate the equilibrium properties of the system. The weighting 
factors p,f of Eq. (5) are all equal to Nf l where Ni is the number of states of the sample. In 
other words, estimates of the equilibrium properties of the system are equal to the average of 
the properties of the states sampled. There are cases where the distribution function , of the 
system is not the optimal choice. Examples are the estimation of properties dominated by the 
contribution of a few rare states (e.g., the estimation of the probability of a molecular system to 
be in a transition state region, or the value at risk of a financial portfolio), or systems that consist 
of several sub-regions and where the sampling gets trapped in one of the sub-regions (e.g., 
simulations of molecular systems, in general, and identification of binding modes of drug-protein 
complexes, in particular). 

The terms "Importance Sampling" and "Umbrella Sampling" refer to sampling algorithms that 
use a sampling distribution function different from the distribution function of the system (MP. 
Allen and DJ. Tildesley, "Computer Simulation of Liquids", Clarendon Press, Oxford (1987)). 
These methods are based on the facts that the distribution function of the system is not always 
the optimal sampling distribution function, and that properties of the unbiased system can be 
estimated from simulations with a sampling distribution different from the distribution function of 
the unbiased system. Importance sampling is a general concept The concept does jiot say 
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which kind of sampling distribution function to use. For complex systems with a large number of 
states, defining appropriate sampling distribution function is very difficult 

US 6,381 ,586 describes a method to get rapid convergence of the estimate of the price of an 
option. The method is based on an appropriate selection of the importance sampling distribution 
function. The method relies on the fact that states that contribute to the estimate are 
concentrated in a sub-region of the state space, and that this sub-region can be identified before 
starting the simulations. If these conditions are fulfilled, convergence of the estimate can be 
increase by identifying the sub-region, and by defining a sampling distribution function that 
biases sampling into the identified sub-region. This situation is quite particular to the problem 
considered, and the same approach is not applicable, for example, for the determination of the 
value of risk of an option portfolio. In many applications it is not possible to identify such a sub- 
region before starting the simulation. Defining a wrong sub-region would cause the simulation to 
be trapped within said sub-region and slow down convergence of estimates. 

US 6,178,384 describes a sampling method to estimate conformational free energies. The 
method is based upon an appropriate definition of the sampling distribution functions. The 
method assumes that regions that contribute to the estimates are approximately harmonic, and 
that the centres of these regions are known. Using these assumptions, a series of Monte Carlo 
simulations is carried out Each of the simulations samples one of the identified regions using a 
sampling distribution function equal to the corresponding harmonic potential. The method is 
difficult to apply, if the system cannot be described as a sum of approximately harmonic wells, if 
there is a large number of harmonic wells that contribute to the estimates, or if the wells cannot 
be identified. For example, molecular systems in explicit water consist of a very large number of 
wells, or the dominant local minima of proteins are difficult to identify without experimental data. 

Defining sampling distribution function to improve the convergence of the estimates requires 
some knowledge of the system. For most systems of interest this knowledge is not available a 
priori, and must be acquired somehow. One possibility is to carry through preliminary studies, for 
example, minimization as is proposed in the patents 6,178,384 B1 and 6,381,586 B1. The 
problem is that with increasing complexity of the system (large number of local minima, many 
dimensions, regions of state space that contribute to the estimates are unknown) the required 
preliminary studies become difficult or impossible. 

Instead of defining the sampling distribution function before starting the simulation, it is in many 
cases easier to identify a desired property of the sampling, and to iteratively adapt the sampling 
distribution function in order to achieve the desired property. This is the essence of what is 
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referred to as "adaptive umbrella sampling" in the literature on molecular simulations (Christian 
Bartels, Michael Schaefer and Martin Karplus "Determination of equilibrium properties of 
biomolecular systems using multidimensional adaptive umbrella sampling", J. Chem. Phys. 111, 
17 (1999) 8048-8067; Bemd A Berg and Thomas Neuhaus "Multicanonical ensemble: A new 
approach to simulate first-order phase transitions", Phys. Rev. Let 68 (1992), 9-12). It can for 
example be assumed that the property of interest is the probability distribution of a selected 
degree of freedom (probability distribution of the distance between two atoms, probability 
distribution of future gains of a financial institution, ...). Estimates of such a distribution require 
that all values of the degree of freedom are sampled. Thus, it makes sense to adapt the 
sampling distribution function such that all values of interest of the degree of freedom are 
sampled with comparable probabilities. A problem with adaptive umbrella sampling is that even 
if the desired umbrella potential is determined, it will in general not be possible to ensure that 
transitions between different important regions of the system occur and that the system does not 
get trapped. 

Summary of the Invention 

The present invention is a method to improve the performance of established sampling 
methods. The new method can be called incremental umbrella sampling. Incremental umbrella 
sampling defines a method for sampling the state space by iteratively generating states xi,t and 
their weighting factors pij by fitting the sampling distribution function pj(x) of the next iteration 
to at least one weighted property of the already sampled states. This means that pj(x) is fitted 
to. the product j&j0(xu), in which the ptj are the weighting factors and 0(pS) is a function 
respectively a property of the states xi,t . 

' In a first step an initial sampling distribution function p(x) is selected. If necessary, other 
starting parameters are selected as well, for example a starting state. The iteration parameter j 
is set equal to a starting value, preferably to 1. The iteration process consists of a second, a 
third, a fitting and a fourth step which are repeated until at least one criterion is fulfilled. In the 
second step Nj states xjj are generated by a numerical sampling algorithm that samples 
states according to the sampling distribution function pj(x) . The third step determines weighting 
factors pij for states xtj generated so far, which means for states xtj with / <>j (the index i 
distinguishes different previous simulations, and the index t distinguishes different states of each 
simulation). The weighting factors pu are determined by using the sampling distribution 
functions p ( (x) used so far, respectively with The number of states » and the number of 



WO 2004/021208 



-6- 



PCT/CH2003/000548 



weighting factors ptj can be incremented with each iteration. For reasons of insufficient 
memory it can be necessary to discard some of the states. In order to have a consistent set of 
weighting factors p it , the weighting factors are calculated in each iteration for all, respectively 

for a set of selected, states. At least for normalized weighting factors p i( the weighting factors 

p ig of iterations with / < j (earlier iterations) will be altered by the newest iteration, respectively 

by the fact that the total number of selected states is increased. The iteration parameter j is 

incremented preferably by 1. the fitting step determines a sampling distribution function pj{x) 

by fitting pj(x) to ptjO{pS) for states xij generated so far. The fitting can be done at all or at 

selected states x tJ . In the fourth step at least one criterion is tested, for example the number of 

iterations. In many cases good results can be found if at least three iterations are done. 
According to the result of the test it will be decided whether to continue the iteration with the 
second or to go to a fifth step in order to perform an analysis. The forth step could as well be 
executed directly after the third step. In this case the iteration would be continued by the fitting 
step and then by the second step. 

By fitting pj(x) in the state space it is possible to use all the information of pj and 0(xu) for 

the states x u generated so far. The fitting step allows to use different fitting strategies. For 

example the fitting can bias the sampling away from areas where intensive sampling has been 
done in the preceding iterations, or the sampling can be directed along local gradients 
respectively towards local minima or maxima of one or several weighted properties. In each of 
the iterations, the sampling distribution function is fitted in a way to improve the overall sampling 
of the state space. For example each fitted sampling distribution function will direct the sampling 
to regions that contribute significantly to the estimates of interest and that were little sampled so 
far. For a given region of the state space the fitted sampling distribution is different if there is a 
large number of states with small weighting factors or if there are a few states with large 
weighting factors. In contrast, with adaptive umbrella sampling methods and other methods that 
are based on estimating some properties of the sampling distribution function with Eq. 5, the two 
situations would result in comparable sampling distribution functions. 

Adaptive umbrella sampling attempts to come up with a single sampling distribution function that 
should ensure sampling of all the important states. Incremental umbrella sampling creates for 
each iteration a sampling distribution function to improve the overall sampling. Since incremental 
umbrella sampling avoids re-sampling of regions that were already extensively sampled and 
directs sampling to regions that were relatively little sampled, incremental umbrella sampling is 
more efficient than adaptive umbrella sampling. 
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The inventive method reduces the number of states that need to be sampled to identify the 
states of interest, or to obtain converged estimates of the state space integrals. Compared to 
simulated annealing (S . Kirkpatrick, C. D. Gelatt Jr., M. P. Vecchi, "Optimization by Simulated 

5 Annealing", Science! 220, 4598, (1983) 671-680) and multicanonical sampling (Bemd A. Berg 
and Thomas Neuhaus "Multicanonical ensemble: A new approach to simulate first-order phase 
transitions", Phys. Rev. Let 68 (1992), 9-12) the invention has several advantages. The 
probability . of the identified optimal state can be estimated. It supports multi-objective 
optimisations. State space integrals can be solved. It reduces the probability that the system is 

o trapped. The invention is general. It can be used with different sampling methods, in particular 
with Monte Carlo sampling, Metropolis Monte Carlo sampling, or dynamic simulations. It can be 
combined with the concepts of simulated annealing and multicanonical sampling. It provides a 
general framework that can be adapted to the system and the observables of interest 

5 A special embodiment uses in at least one iteration a numerical sampling algorithm which 
generates correlated states xjj , wherein the sampling distribution function pj(x) preferably has 
a maximum that biases sampling into a sub-region of the state space, wherein the sampling 
distribution function pj(x) is preferably fitted such that the maximum is in a region where the 
product pijOipj) has a maximum, and wherein the sampling preferably starts from a state 

:o close to the maximum of the sampling distribution function pj(x) . 

The inventive method is preferably commercialized in the form of a computer software product 
on a computer-readable medium in which program instructions are stored, which instructions, 
when read by a computer, enable the computer to apply the inventive method. 

:5 

Brief Description of the Drawings 

Fig. 1. is a flow chart of an iterative umbrella sampling simulation. 

Fig. 2. shows a distribution function p(x) (asterix), an observable 0(x) (circles) and a product 
o &)p(x) (rectangles) of a one-dimensional system in the form of probabilities (p) versus the 
dimension of the system (x). 

Fig. 3. shows distribution functions pi(x) (circles), pi{x) (rectangles) and f*{x)/2+ pi{x)/2 
(asterixes) of the one-dimensional system (p versus x). 

Fig. 4. shows the weighting factors /5m (circles) and the product putd{xi % i) (asterix) versus the 
•5 states xt\t . 
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Fig. 5. is an example of a conformation of the two-dimensional HP protein-folding model. 
Fig. 6. shows estimates of the probability (p) of states with a given energy (E) for the HP protein- 
folding model. 

Fig. 7. shows estimates of the probability of states with energy -92, -90, -88, -86, -84 and -82 as 
a function of the number of states that were sampled for the HP protein-folding model. 
Fig. 8. shows energies (E) of states sampled in test run B as a function of the number of states 
(n) that were sampled. 

Fig. 9. shows estimates of the probability of states (In p) with given radius (R) for runs A, B, C, D 
and E. 

Fig. 10. shows energies (E) of states sampled in the simulated annealing run B as a function of 
the number of states (n) that were sampled. 

Fig. 11. shows energies (E) of states sampled in the multicanonical (Adumb) run B as a function 
of the number of states (n) that were sampled. 

Fig. 12. shows estimates of the probability (p) of states with a given energy (E). Estimates 
obtained from the five multicanonical (Adumb) runs A, B, C, D, E and the average of the test 
runs from Fig. 6 Ref. 

Detailed Description of the preferred Embodiments 

Exploration of the state space proceeds in a series of steps shown in Fig. 1. In a first step 1 an 
initial sampling distribution function^>c>i(x:) is selected. If necessary other starting parameters are 
selected as well, for example a starting state. This is often an interactive step. Computational 
tools may be used. The iteration parameter j is set equal to 1 . The iteration consists of a 
second 2, a third 3, a fitting F and a fourth step 4 which are repeated until at least one criterion is 
fulfilled. In the second step Nj states xu with f=l,2,.„,iV) are generated by a numerical 
sampling algorithm. In the third step 3 the weighting factors p it are deduced from the sampling 
distribution functions p t (x) for states xu generated so far, which means for states xti with / 
<, j . The third step 3 increments j by 1. The fitting step F determines the sampling distribution 
function pj(x) for the next iteration by fitting pj(x) to puOipcu) for states xu generated so far. 
In the fourth step 4 at least one criterion is tested. For example the criterion can be the 
convergence of the simulations or reaching a given number of iterations. The criterion may 
involve some interactive steps. According to the result N or Y of the fourth step 4 it will be 
decided whether to continue (N) the iteration with the second step 2 or to go (Y) to a fifth step 5 
in order to perform the final analysis. The simulations, respectively the states xj,t with the 
weighting factors p Ut , are analysed to obtain the estimates of interest An analysis based on 
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Eq. (5) is convenient, as it provides a general and consistent way to obtain estimates of 
equilibrium properties. 

The initial conditions should be selected such that the simulation reaches rapidly equilibrium. 
This is necessary for the sample to represent the sampling distribution function. If a Metropolis 
Monte Carlo algorithm or a molecular dynamics simulation is used, the starting state and the 
initial sampling distribution function should be selected such that the starting state has a high 
probability of being sampled with the selected sampling distribution function. If a Monte Cario 
algorithm is used that produces sates independently of any starting state, selection of a starting 
state is evidently not required. 

A simulation is carried out with the selected sampling distribution function pj(x) . The simulation 
must produce a sample consisting of JV) states xjj distributed according to the sampling 
distribution function pj{x) . 

Different simulation algorithms may be used, for example Monte Carlo, Metropolis Monte Cario, 
Molecular Dynamics algorithms. The choice will depend on the system to be studied and the 
selected sampling distribution function. For reference in later sections, some possibilities are 
listed here 

a) Monte Cario Method: States are generated independently of each other by 
transforming random numbers generated by a random number generator. The resulting 
states are un-correlated. 

b) Metropolis Monte Carlo, Version 1: A trial state xtnai is generated from an 
existing state xu such that the probability of generating xtnai from xjj is equal to the 
probability of generating xjj from xtnai . The trial state is accepted with 

probability 1, if pfaiafepfejj) (6A) 

probability pi{xtriai)f pi(xj t t) if p{xtrtai)<pi{xjj) (6B) 

If xtnai is accepted, it is taken as the new state xit+i=xtriai , otherwise xjj is taken as the 

new state xu+\=xjj . 

c) Metropolis Monte Carlo, Version 2: Using a Monte Carlo algorithm that 

produces states distributed according to fhim P i^x) t a trial state xtnai is generated 

v 

independently of any other state. The trial state is accepted with 

probability 1 , if v> , (7A) 

J Pslmple{Xtrial) p,imple\X] % t) 
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probabjlity P'Mpwixj,) jf pfa*) pM (7B) 

r Pslmple\Xtrlal)pi\Xjj) p»mple\Xlrial ) pilmpU\X),t) 

If .ttrt./ is accepted, it is taken as the new state xjj+^xiriat , otherwise xj,t is taken as the new 
state Xj.M=xj.t . 

New sampling methods are still being developed. US 5,740,072 describes a sampling method 
5 that combines the Metropolis Monte Carlo method with stochastic dynamics. 

Optimal weighting factors 



Optimal weighting factors can be obtained with Equations (8-11). Equations (8) and (9) express 
io the weighting factors pff in function of the normalisation constants /{# , and vice versa. 

Self-consistent estimates of the weighting factors pff and normalization constants f k self can 
be obtained by applying equations (8, 9) repeatedly until the weighting factors have converged. 
15 The normalised weighting factors pv are then given by 

In equations (8) to (10), Nj and AT* are the number of states produced by simulations; and k, 
respectively, Nsim is the number of simulations, and the bias functions 

c*{x)*:p k {x)lrtx) < 11 > 
>o are proportional to the ratio of the sampling distribution function of the simulation divided by the 
distribution function of the system. It helps with the understanding of what follows to note that the 
weighting factors ptj are proportional to pfa,t)f , if only a single simulation is analysed. If 
several simulations are analysed, the weighting factors can be understood as the ratio of the 
distribution function of the system divided by an average effective sampling distribution function 
15 that represents the sampling carried out in all of the simulations. Thus, the weighting factors are 
large for states, for which the distribution function of the system is large, or for which the average 
effective sampling distribution function is small. A necessary condition for any method to 
estimate the weighting factors to be useful in the context of the present invention is that 
additional sampling in the region around a selected state must reduce the weighting factor of the 
30 selected state. 
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For particular cases, other equations exist to determine the set of weighting factors. If the 
sampling distribution functions differ only along one or a few selected degrees of freedom, the 
Weighted Histogram Analysis Method of A.M. Ferrenberg and R.H. Swendsen, "Optimized 
Monte Carlo Data Analysis", Phys. Rev. Lett. 63, (1 989) 1 1 95- may be used. If the normalisation 
constants of the sampling distribution functions are known and need not to be estimated with 
Eq. (9), the formalisms described and referenced in Art B. Owen and Yi Zhou, "Safe and 
effective importance sampling", Journal of the American Statistical Association 95 (2000) 135- 
may be used. 

One-dimensional system to illustrate the inventive method 

A fitting procedure is used to define a sampling distribution function that directs sampling to 
regions that contribute significantly to the estimates of interest and that were little sampled so 
far. One of the advantages of the proposed method is that it can be applied to high-dimensional 
systems with sampling distribution functions that vary along many dimensions. The inventive 
method is illustrated using a one-dimensional system, e.g., two atoms separated by a distance 
x. The distribution function of the system of the example is equal to a normal distribution with 

mean 0 and standard deviation 1, i.e., p(x)=-j±=e "* ^ . The observable assumed to be of 

interest is set equal to a normal distribution with mean 2 and standard deviation 1, i.e., 
Q(x)=^=e" (jc " 2 ^ 2 . Further, it is assumed that two simulations have been carried out, and that 

the sampling distribution function for the third simulation should be determined, i.e., j=3 . 

Fig. 2 illustrates the distribution function of the system and the observable; the asterix show the 
distribution function p(x) 9 the circles the observable Ch(x) , and the rectangles show the product 
p(x)Oi(x) used in Eq. (1). The product p(xfii{x) is maximal at x=l. Thus, to calculate the 
average (a) of the property Ch(x) 9 states around x=l are most important. Fig. 3 shows the 
sampling distribution functions pi(x) as circles, pz(x) as rectangles, and the average 
pi(x)/2+pz(x)/2 as asterixes, respectively. The sampling distribution function pi(x) was chosen 
to be equal to 1 for x between -0.5 and 0.5, the sampling distribution function pz(x) was chosen 
to be equal to 0.5 for x between 0 and 2. In each of the two simulations, 30 states were 
sampled. 

Fig. 4 plots for each state xu the weighting factors pu and the product ffoOifo) as circles 
and asterixes, respectively. The sampling distribution functions are normalised. Thus, the 
normalisation constants need not to be estimated with Eq. (9), and the (un-normalised) 
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weighting factors are directly given by 

Tf-psffl&fi' (12> 

They are equal to the distribution function of the system divided by the average of the sampling 
distribution functions of the two previous simulations. 

The estimate of the average (Oi) is dominated by states for which the product pjfh{xij\ is 
large (Eq. 5). From Eq. (12) and Fig. 4, it is seen that the product pi^Gipa^ is large for states 
with ffa t t\Oi(xiJ\ large and p{xiS)/2+pi(xi % t)l2 small. f{xJ\Ch(xiJ\ is large for states that 
contribute significantly to the average (Ch) of the property Oi(x) (Eq. (1) and Fig. 2). 
pi(xij)/2+pi(x!jy2 is small for states that are in regions that were not much sampled so far (Fig. 
3). Thus, to sample states that contribute to the estimate of the average (Ol) of the property 
Oi(x) , and to sample regions that were little sampled so far, sampling should be directed 
towards regions around states that correspond to a large value of the product pv|Ol(*v)| . This 
can be done by fitting the sampling distribution function /a(x) for the next simulation to the 
product pi)fh{xi4\ evaluated for the states obtained from the first two simulations. The fit is 
done such that the resulting sampling distribution function is large for at least some states with a 
large value of the product pi t t\Oi(xi t t] . 

Formulating this more generally: the sampling distribution function Pj(x) for the next simulation 
is fitted to the product pjOfet) of the states of all previous simulations i=l,2 r .^/-l times a 
function 0(x). The fit is done such that the resulting sampling distribution function is large for at 
least some states with a large product puOfa} , and tends to be small for states with a small 
product puOfa). The function 0(x) is defined such that it is large for the states of interest, 
e.g., it is set equal to (X*)=|Cfc(*)| . 

Different strategies can be used to fit the sampling distribution function for the next iteration to 
the product p^Ofa] shown as asterixes in Fig. 4. A simple strategy is to identify the state with 

the largest value of the product pi\Oi{xtJ\ , and to define the sampling distribution function for 
the next iteration such that it is large at the selected state. In the example of Fig. 4, the state with 
the largest value is the 16 th state of the second simulation and is located at 1.07, i.e., **.i6=1.07 . 
Thus, a possible sampling distribution function for the next iteration of the example might be 
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uniform sampling around the state at 1.07 in the range between 0.57 and 1.57. An alternative 
strategy would be to identify states with a small value of the product pi\Ch(xij\ and to make the 
sampling distribution function for the next iteration small for these states. E.g., the sampling 
distribution function for the next iteration of the example could be set to a uniform distribution 
between -0.5 and 2 but excluding the region from -1.37 to 0.63 that lies around the state at - 
0.37 that has the smallest value of the product pjfk{xij\ . A third strategy would be to express 
the logarithm of the sampling distribution function as a linear function of some parameters and to 
use a linear least square fit to determine the parameters. These are only some of the possible 
fitting strategies. They are described in a more formal way below. 

Formal description of a fitted sampling distribution function 

The sampling distribution function pj(x) for the next iteration is defined based on the states 
j^and their weighting factors pu from the previous simulations M^.-J-l. The sampling 
distribution function pj(x) can be expressed as a function h(yj,x) of a vector vj of parameters, 
and the normalisation constant fj , i.e., 

pj{x)=pj(yj,x)=fj h(yj,x) . (13) 
For some of the existing simulations algorithm, e.g., molecular dynamics or Metropolis Monte 
Carlo, the normalisation constant fj needs not to be known, and needs not to be determined. 
For the other simulation algorithms, it can be determined by applying the constraint that the 
integral of the sampling distribution function over all states must be normalised to one. 

The parameters of the sampling distribution function are determined by fitting the sampling 
distribution function p/x) to the weighting factors pu from previous simulations ze{l,2,...,y-l} 
and *e{l,2,...,iV/}. The fit maximises an objective function G . The objective function G is 
defined as a function of local comparisons of the sampling distribution function pjiyjyx) with the 
product jxjOfeu) of the weighting factors pu times a function 0(x). The objective function G 
must be large for sampling distribution functions that are large for at least some states with a 
large product fx,tO(xu) , and small for the majority of states with a small product p t tO{xu) . 
Different procedures can be used to do the fitting; examples are given below.. 

Objective functions G are defined as functions of local comparisons. For the purpose of the 
local comparison, states 5Sv from the previous simulations f=l,2,...,/-l are associated to the 
weighting factors p,t from previous simulations, e.g., %j=xu . The sampling distribution function 
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p/yjjQj) at 55./ is compared to the corresponding product pijO{xij) using a local comparison 
function 

du^(pj(yj,xs)&j,pijO(xu)) , (14) 
e.g., rf(0/(vs,5^ , and the objective function is defined 

as a function of the local comparisons, e.g., {^2*72? ditt ' 

The form of the sampling distribution function, the set of parameters, the objective function G , 
the function 0(x), and the fitting procedure may be adapted to the system and properties of 
interest, as well as to the available simulation algorithms. Examples are given below. 

Examples of sampling distribution functions with hanmonic constraints 

Examples of sampling distribution functions can be found by applying harmonic constraints to 
the distribution function of the system, 

Pj{x)= fjGxJ^kcoitstr{x-Xj t constry )fi{x) . (15) 

In Eq. (15), kco^tr is a constant that determines the strength of the harmonic constraint, fj is 
the scaling factor to normalise pj(x), and xj.comtr is the centre of the constraint The parameters 
that can be fitted are the centre and the strength of the constraint, vj={kconstr 9 x J \constr} . 
For systems consisting of many local minima separated by barriers, application of a 
multicanonical sampling distribution function can be advantageous, 

Pj(x)= fjGX^kconstr(x-Xj t constry)pmc(x) . (16) 

In Eq. (16), p mc (x) denotes the multicanonical distribution function; the remaining parameters 
are identical to those of Eq. (15). In particular, the parameters that can be fitted are the centre 
and the strength of the constraint, vj={kconstr,xj,constr} . 

The multicanonical distribution function pmc(x) (Bemd A. Berg and Thomas Neuhags 
"Multicanonical ensemble: A new approach to simulate first-order phase transitions", Phys. Rev. 
Let 68 (1992); Christian Bartels, Michael Sqhaefer and Martin Karplus "Determination of 
equilibrium properties of biomolecular systems using multidimensional adaptive umbrella 
sampling", J. Chem. Phys. 111, 17 (1999) 8048-8067) is determined such that different values 
of ln/?(jc) are sampled with comparable probability. The multicanonical distribution function can 
be written as 

/^X^JCy/mrOn/^)) , (17) 

where f mc denotes an estimate of the probability distribution of hip(x) . The estimate fmc can 
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be obtained from a histogram of lnp(x) that can be obtained from the sampled states^ and 
their weighting factors pu using Eq. (5) . 

A sum of harmonic constraints applied to the distribution function of the system results in the 
distribution function 

pj{x)= fjeXT^"'^ kjJ t constr(x-XjJ,constr j)rfx). (18) 

In Eq. (18), the constants kjj t constr determine the strengths of the harmonic constraints, fj is the 
scaling factor to normalise pj(x) f xj^constr are the centres of the constraints, and N^con** is the 
number of harmonic constraints. The parameters that can be fitted are the number of harmonic 
constraints, their centres, and their strengths, 

Vj~{kj t i .constrykjXconstr^ . ^Xj,l t constr 9 Xj2 t constr>. . Nj t constr } . 

Depending on the sampling distribution function pj(x) 9 particular sampling methods may be 
more or less useful. For example, the sampling distribution functions defined by Eqs. (15) or (16) 
with a positive value for the constant kcomtr constrain sampling to a sub-region of the state 
space and are particular useful together with sampling methods such as the Metropolis Monte 
Carlo Version 1 or stochastic dynamics that generate correlated states. In contrast, if the 
sampling distribution function of Eq. (18) is used, convergence of a Metropolis Monte Carlo 
Version 1 algorithm may be slow, and Metropolis Monte Carlo Version 2 or Monte Carlo 
algorithms may be preferable. 

Examples of the function 0(x) 

Examples of the function 0(x) can be defined such that p(x)o{x)\s large for states that should 
be sampled in one of the simulations. The states that should be sampled are primarily those for 
which one of the N P ro P properties Q^a(x}(h(x\.J)»„(x)} (see below and Eq. 1) is large. 
The function 0(x) is set equal to a function of the set 0 of properties, e.g., 

Ofa max l^M (19) 
or 0(x)=^ l^fe^ with the normalisation constant sk for each property k defined as 

Similar as in umbrella sampling, three different types of properties can be distinguished that 
identity important states and should be included in the set 0 of properties. First, states may be 
important since they occur with large probability in the unbiased system (p(x) is large): In this 
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case, a property Ok(x)=l should be added to the set © of properties. Second, states may be 



the final analysis. Such a property should be added to the set 0 of properties. Third, there may 
exist states that are not directly of interest but whose sampling is important to obtain converged 
estimates of other properties. For example, in molecular systems where one is mainly interested 
in low energy conformations, sampling of high-energy conformations is advantageous to ensure 
transitions between different regions of low energy. In this case, a function should be added to 
the set © that is large for states that need to be sampled to increase convergence of the 
estimates. Examples are given below. 

Fitting Procedure 1: Fit to state with large weighting factor 

Procedure 1 fits a sampling distribution function such that it is particular large for a selected state 
with a large product puOipaJ). The description that follows assumes that the sampling 
distribution function of Eq. (16) is used, that the centre xuons* of the constraint is determined by 
the fit, and that the value of the constant kconstr is set to a positive value that is determined in a 
series of preliminary test runs. The procedure may be adapted, e.g., for other sampling 
distribution functions. 

Based on the value of the product p.<0(*?.') . the state x*a from the already sampled states is 
identified that has a comparatively large weighting factor. E.g., the criterion to identify the state 
xsei can be that it corresponds to the maximum 



In Eq. (20), the maximum is taken over the previous simulations /, and over the states t 
generated in each of the simulations; N*im= y-1 is the number of simulations earned out so far. 
The centre of the constraint is then set equal to the selected state, i.e., xj tC onstr=xsei . The 
selected state xsei may also be used as the starting state for the simulation of the next iteration. 

The objective function that is maximised by procedure 1 is equal to 



important since they contribute significantly to an average (Ot) of a property Ok(x) of interest in 



,max . max JSuOfpaj). 



(20) 




(21) 



(22) 
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Fitting Procedure 2: Fit to states with small weighting factors. 

Procedure 2 fits a sampling distribution function to a state with a particular small weighting 
factor. The description assumes again that the sampling distribution function of Eq. (16) is used 
and that the centre xj.a>nstr of the constraint is determined by the fit The value of the constant 
kcomtr is set to a negative value that may be determined in a series of preliminary test runs. The 
procedure may be adapted; e.g., for other sampling distribution functions with a defined centre. 
A state xsd is selected with a small product JktOfat) , e.g., the one that corresponds to the 
minimum 

t «&~),A}fr°(*<y (23) 

The centre of the constraint xj^omtr is set equal to the selected state x$ei . The objective function 
that is maximised by the procedure can be written as 

G=- ,min . inin ///,r , (24) 

where the local comparisons are defined by Eq. (22). 

More generally, a sampling distribution function can be fitted to a set of states with particular 
small weighting factors, e.g., the centres xjxeomtr may be determined of the sampling distribution 
function of Eq. (18) with kjuonstr set to negative values. For this purpose Nj^mstr states Xseu 
with j=1,2,. . .,Nj,constr may be selected that correspond to the minima 

,- 6{ ,»^,e{& (25> 
This may be done by selecting first the state x*eu that corresponds to the minimum of Eq. (25) 
with 5=1 , then select the state Xseia with 5=2 , and so on until s=Nj t constr . The centres, xjj^mr , 
of the constraint can then be set equal to the selected states, xseij . This procedure maximises 
the objective function of Eq. (24) with the local comparisons defined as 

dij^(j?j(vj 9 xj)%j 9 pijO(xij))= 

Fitting Procedure 3: Linear least square fit 

Procedure 3 performs a linear least square optimisation to fit the vector vj of parameters of the 
sampling distribution function to the product psJf{pS\ . The description assumes that the 
sampling distribution function of Eq. (18) is used. The number of harmonic constraints, Nuonstr , 
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and the centres, xjuonur, of the harmonic constraints are determined before starting the fit This 
can be done, for example, using a clustering algorithm. The strengths kjAcomtr of the constraints 
are obtained by a linear least square fit of the logarithm of the un-normalised sampling 
distribution function to the logarithm of the product pi t \o{xiJ\ of the weighting factors ptj times 
the absolute value of the function 0(x). The objective function can be set equal to 



The procedure may be adapted, e.g., for other sampling distribution functions. 

Fitting Procedure 4: Fit to pairs of states with large and small weighting factors, respectively. 

Procedure 4 fits the sampling distribution function to pairs of states x S ei\ and xsen with the two 
states being close together in the state space, and with p*en being relatively small and psea 
being relatively large. Such a pair can be identified, for example, using the criteria of Eq. (20) to 
select xsen , and then select the state closest to it as xseix . An alternative is to select states 
generated in subsequent time steps with a Metropolis Monte Carlo Version 1 , or a stochastic 
dynamics algorithm. The selected states are then used to define a bias away from xsen into the 
region around xsen , e.g., use the sampling distribution function of Eq. (16) with 
xuconstr = Xsen+cfxsen-xsdi) , where c is a constant 

The objective function maximised by this procedure is the same as the one maximised by 
procedure 1 , except that the state at which to evaluate the sampling distribution function should 
be defined as xj^xu+cfaj-xiAciose) , where jcucf«» denotes the state from the sampled states 

that is closest to xu . 

Sampling of conformations of molecules 

Sampling molecular conformations is one of the main applications of molecular dynamics 
simulations. Examples are peptide and protein folding simulations, in which #ie conformational 
space of a peptide or protein is investigated with the aim of identifying conformations that are 
stable under native conditions, or drug binding studies with the aim of identifying and 
characterising binding modes of drug candidates to their host proteins. As an example, the 
following description describes how to apply the patent to the protein folding problem. 



(27) 



with local comparisons defined by 



(28) 
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Step 1. Perform a short molecular dynamics simulated annealing (S . Kirkpatrick, C. D. Gelatt 
Jr., M. P. Vecchi, "Optimization by Simulated Annealing", Science, 220, 4598, (1983) 671-680.) 
run starting from an extended conformation. The conformation, Xm^i, resulting from the 
annealing run is used as the starting state of the simulation when executing step 2 for the first 
time. In the first simulation, the system is restrained to the region around x«w , i.e., the 
sampling distribution function of Eq. (15) is used with xx.cn^xmneai . The parameter ken,* that 
determines the strength of the harmonic constraint is set to a positive value that can be set 
based on a set of preliminary test runs. 

Step 2. Carry through a molecular dynamics simulation starting with the selected starting state 
(conformer), and the selected sampling distribution function pj{x) . 
Step 3.1. The weighting factors pi.t are determined with Eqs. (8-10). 

Step 3.2. Fitting Procedure 1 is used to fit the sampling distribution function of Eq. (16). A single 
property is used to define the function C{x) with Eq. (19) and to select a state with Eq. (20). The 
property is set equal to the inverse of the probability distribution of the energies of the system. 
The probability distribution of the energies is estimated using Eq. (5) based on the energies of 
the sampled conformations xu and their weighting factors p.t. The property is large for 
conformations with energies that are unlikely in the unbiased system, and it is small for 
conformations with energies that are likely in the unbiased system. Using this property in Eq. 
(19) ensures that states are selected regardless of their energy, provided they are in a region 
that was relatively little sampled in previous simulations. Depending on the system, additional 
properties could be included in the set 0 of properties used to define the function 0(x). For 
example, for drug binding studies, the additional property might be set equal to the inverse of the 
probability distribution of the position of the drug candidate relative to the protein. This would 
ensure that different possible positions of the drug relative to the protein are sampled. 
Step 4. Make 200 iterations of steps 2 to 4. 

Step 5. Determine the weighting factors job. Use the weighting factors to estimate properties, 
e.g., the free energy as a function of the temperature, or the probability distribution of the radius 
of gyration at a selected temperature. 

Two-dimensional HP nrotein-foldinq model 

This section applies the invention to study the two-dimensional HP protein-folding model, and 
compares it to simulated annealing and to multicanonical sampling. The conformational space is 
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investigated of a chain consisting of 100 hydrophobic H and polar P monomers. The chain 
moves on a two-dimensional grid. The chain is self-avoiding, i.e., no two monomers may occupy 
the same grid point If two H monomers are in contact, they interact favourably giving an energy 
contribution of -2. There are .no other interactions. The probability distribution of the system is 
set equal to fa)=exz(-E(x)/2) where E(x) is the energy function. The probability distribution 
is not normalised. A Metropolis Monte Carlo algorithm is used to sample conformations of the 
chain. The Metropolis Monte Carlo algorithm uses a set of simple moves to change the 
conformation of the chain locally, one or two monomers per move. The same set of moves is 
used in all the simulations, and no attempts were made to optimise the set of Monte Carlo 
moves. Fig. 5 shows a conformation. Hydrophobic H monomers that interact favourably if in 
contact are drawn as empty rectangles; polar P monomers are drawn as filled circles. 

To study the conformational space using simple enumeration, some 1.7*10 47 conformations 
would have to be generated, which is beyond the capabilities of computers. The examples 
presented in this study change the conformation 3*10® times. The required calculations take a 
few hours on a personal computer. The properties that are assumed to be of interest and that 
are subsequently analysed are low energy states, states with small radii, the probability of states 
as function of the energy and as a function of the radius. The first two properties illustrate global 
optimisations; the latter two properties illustrate state space integrals. 

The procedure used to study the folding of the grid model is similar to the one described above 
for the protein folding simulations. In particular, a sampling distribution function similar to the 
distribution function of Eq. (16) is used that applies a local constraint, and the sampling 
distribution function is fitted using fitting procedure 1. In more detail: the constraining potential is 
defined based on the radius of gyration and on a comparison of the monomer contacts present 
in the current conformation x and the conformation selected by the fitting process xsei . I.e., the 
sampling distribution function is set equal to 

pj(x)=fjext{-fo(nformed(x)±mn^ . (29) 

In Eq. (29), nformed(x) is equal to the number of monomer contacts formed in the current 
conformation xthat are not present in x** t itbrokJjc) is equal to the number of monomer 
contacts of xsa that are no more present in the current conformation x . Rgyr(x) is the radius of 
gyration of conformation x , and Rref is set equal to the radius of gyration of x*ei . Since a 
Metropolis Monte Carlo algorithm (version 1) is used to generate the conformations, the 
normalisation constant fj in Eq. (29) needs not to be known. The constants h and ki are set 
equal to 0.003 and 0.5, respectively. The constants were determined in a series of test runs. 
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For the simulations (step 2), a simple Metropolis Monte Carlo algorithm (version 1) is used. The 
Monte Carlo algorithm generates new conformations by randomly changing the position of a 
single monomer or two adjacent monomers at a time. New conformations are accepted or 
rejected using the criterion of Eq. (6). A single run of the simulation algorithm in step 2, consists 
5 of 3'000'000 Monte Carlo moves. The first 500*000 moves are considered as equilibration 
phase, and are not used. Then, 50 conformations are taken every 50*000 moves. This results in 
50 conformations xu for each of the simulations j. 100 simulations are carried out. Thus, a total 
of 3*1 0 B states are generated, of which 5*000 are used for the analysis, e.g., to fit the sampling 
distribution functions in step 3. 

o 

Two properties are used to define the function 0(x) with Eq. (19) and to fit the sampling 
distribution function using procedure 1. The first property is the inverse of the probability 
distribution of the energies of the system. The probability distribution of the energies is estimated 
using Eq. (5) based on the energies of the sampled conformations xu and their weighting 
5 factors pu. The property is large for conformations with energies that are unlikely in the 
unbiased system described by p(x), and it is small for conformations with energies that are 
likely in the unbiased system. The property ensures that states are selected regardless of their 
energy, if they are in a region that was little sampled in previous simulations. 

0 The second property used in Eq. (19) is the inverse of the probability distribution of the radius of 
gyration of the system. This property ensures that states are selected regardless of their radius 
of gyration, if they are in a region that was little sampled in previous simulations. The selection of 
the properties to use for the fitting procedure is compatible with the set of properties that are 
assumed to be of interest (low energy states, states with small radii, the probability of states as 

5 function of the energy and as a function of the radius; see above). 

Results from five independent test runs are analysed. Fig. 6 plots the probability of states as a 
function of the energy. The agreement between the results of the different runs is good. Some 
differences exist for the estimates of the probabilities of states with low energies. The deferences 
o in the estimates suggest that there exist different kinds of low energy states, and that the runs 
differ in the set of low energy states sampled. Low energy conformations with energies around - 
90 are rare, i.e., they are estimated to occur with a probability of about 10* 45 

Fig. 7 illustrates the convergence of the estimates of the probability of low energy states from 
5 Fig. 6. Shown are the estimates of the probability as a function of the progress of the simulation. 
Low energy conformations are rare, and finding them is difficult They are found only after some 
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initial sampling that explores the state space. Once, a first state with given energy has been 
found, the estimate of the probability converges rapidly, and does remain stable over the rest of 
the run. 

5 Fig. 8 plots the energy of the states from Run B as a function of the progress of the run. The run 
starts by sampling high energy states. As the run progresses, it finds states of increasingly low 
energy. The system gets never trapped. After having found a state of low energy, the run 
continues by sampling new regions of the state space, and eventually finds a new low energy 
state. 

10 

Fig. 9 shows estimates of the radius distribution of states. The estimates from the different runs 
are all smooth functions and agree quantitatively. 

Comparison with the state of the art 

15 

In summary, the five runs are successful in identifying low energy and compact states, and give 
quantitative estimates of the probability of the energy of states and of the radii of states. The five 
test runs are compared to five simulated annealing runs and to five runs of multicanonical 
sampling. The same Metropolis Monte Carlo algorithm is used in all the runs, and, as in the test 

20 runs, the conformations are updated 3*1 0 8 times. In the simulated annealing run (S . Kirkpatrick, 
C. D. Gelatt Jr., M. P. Vecchi, "Optimization by Simulated Annealing", Science, 220, 4598, 
(1983) 671-680.), the temperature is lowered linearly in 5000 steps from 1.5 to 0.5. In the 
multicanonical runs (Bemd A. Berg and Thomas Neuhaus "Multicanonical ensemble: A new 
approach to simulate first-order phase transitions", Phys. Rev. Let 68 (1992), 9-12.), the 

15 umbrella potential is updated 1 00 times. 

Table 1 compares the runs with respect to their performance for the optimisation problem "find 
low energy states". The lowest energy found in different runs varies between --80, and -94. 
From all the runs, the test runs with the present invention identify the best state (energy of -94). 
30 In average, the test runs find states with an energy of --88. This compares favourably to the 
average of -86.4 from the multicanonical runs, and to the average of -83.6 from the simulated 
annealing runs. 



Min Score 


Run A 


Run B 


Run C 


Run D 


Run E 


Adumb 


-86 


-86 


-86 


-86 


-88 


Annealing 


-86 


-88 


-84 


-80 


-80 


Test 


-86 


-90 


-90 


-80 


-94 
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Table 1 : Lowest energies found in each of the runs. Data is shown for the five test 
runs (Test), the five multicanonical or adaptive umbrella sampling runs (Adumb), 
and the five simulated annealing runs (Annealing). 

Runs with the present invention do not get trapped (see Fig. 8). This is different from the 
situation with simulated annealing and multicanonical simulations. In simulated annealing runs 
(see Fig. 10), the simulations start at a high temperature at which transitions between different 
regions occur frequently. As the temperature is lowered, sampling is increasingly restricted to a 
region around a low energy state, if different regions with comparable low energies exist - as in 
the present system .- several simulated annealing runs are required to identify the different 
regions. 

Multicanonical simulations of high dimensional system have the tendency to get trapped (Fig. 
11) close to the first low energy state found. Although multicanonical simulations should 
guarantee that low and high energy states are sampled with comparable probability, they cannot 
guarantee that transitions between the different regions occur frequently. 

Fig. 12 shows estimates of the probability of states with given energy from the multicanonical 
runs and compares them to the corresponding estimates from the test runs (Fig. 6). The 
average of the estimates from the test runs is shown as continuos line in Fig. 12. The estimates 
from all the runs agree in general. The estimates for high-energy states agree well. Some 
differences are evident at low energies. At low energies, the probabilities estimated from the test 
runs are larger than those from the multicanonical runs. This indicates that the test runs find 
more low energy states. 

The simulated annealing runs are irrelevant in such a comparison of equilibrium properties of 
the system, since no method has so far been described to derive converged state space 
integrals from simulated annealing runs. 

Value at Risk of a Financial Portfolio 

The value at risk of a financial portfolio gives an estimate of the future value of a financial 
portfolio in adverse market situations. It is equal to a percentile of the probability distribution of 
the future value of the portfolio. E.g., the 99% value at risk is equal to the future value of the 
portfolio that separates the worst 1% from the best 99% of the possible future values of the 
portfolio. Efficient estimation of the value at risk is challenging since it requires sampling of a few 
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rare state, i.e., particular adverse market situations, and the location of the rare states is, in 
general, not known before starting the simulations and must be found. 

Step 1. Select the distribution function of the (unbiased) system as the initial sampling 
distribution function, i.e., p(x)^/tx). A Monte Carlo algorithm is used for the simulations. The 
Monte Carlo algorithm produces states independently of any starting state. Thus, no starting 
state is needed. A Monte Carlo algorithm must exist to produce un-correlated states distributed 
according to /a(x). This restricts the set of sampling distribution functions that can be used. 
Step 2. Cany through a Monte Carlo simulation using the selected sampling distribution function 

> pj(x). 

Step 3.1. The weighting factors can be determined with Eqs. (8) and (10). The normalisation 
constants f k sd/ are known and need not to be determined with Eq. (9) 

Step 3.2. Fitting Procedure 1 is used to fit the sampling distribution function of Eq. (15). A single 

> property is used to define the function 0(x) with Eq. (1 9) and to select a state with Eq. (20). The 
property is defined as follows. First, determine a limit identifying large losses, e.g., losses that 
occur with a probability of less than 10%. Then, set the property equal to a function that is 0 for 
small losses and 1 for large losses. 

) Step 4. Make 200 iterations of steps 2 to 4. 

Step 5. Estimate the weighting factors Use the weighting factors to determine the 
probability distribution of the value of the portfolio, and determine the value at risk. 
An alternative would be to use the sampling distribution function defined by Eq. (18), the 
Metropolis Monte Carlo Version 2 algorithm, and the fitting procedures 2 or 3. 
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Claims 

A Method for sampling a state space by iteratively generating states xtj and their 
weighting factors /5>v , wherein the index / is the iteration parameter and the index t 
distinguishes different states xu generated by an iteration /, the method is consisting of 
a first step for selecting an initial sampling distribution function pi(x), 
a fifth step for performing an analysis 
and an iteration procedure including 

a second step for generating Nj states xu by a numerical sampling algorithm and 
a fourth step for testing at least one criterion to decide whether to continue the iteration 
procedure or to stop the iteration procedure and to go to a fifth step in order to perform the 
analysis using the simulated data, 

characterized in that the iteration procedure further includes 

a third step determining weighting factors pu for states generated so far by using 

sampling distribution functions p ( {x) determined so far and 

a fitting step for determining a sampling distribution function p } (x) for the next iteration by 
fitting pj(x) to jijOfa) for states xu generated so far, wherein 0(xu) is a function, 
respectively a property, of the states x»j . 

Method as claimed in claim 1, wherein the sampling distribution function pj(x) of at least 
one iteration is fitted such that it maximises an objective function preferably defined as a 
function of local comparisons between the sampling distribution function and the product 
jxtO{xi t t) . 

Method as claimed in claim 1 or 2, wherein the sampling distribution function pj(x) of at 
least one iteration is fitted such that the sampling distribution function pj{x) is large for at 
least one state x*,t with a large product puO{xu) , and tends to be small for states with a 
small product puO{xu) . 

Method as claimed in one of claims 1 to 3, wherein the sampling distribution function 
pj(x) of at least one iteration is a function with at least one constraint, preferably having 
an extreme value in the region of at least one selected state xj^tr and having vanishing 
values at a distance from the at least one selected state xj^tr , wherein the constraint is 
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preferably a harmonic constraint with at least one constant keomtr defining the strength of 
the constraint at the selected state xj,constr . 

Method as claimed in claim 4, wherein the sampling distribution function pj(x) of at least 
one iteration is the distribution function of the system with constraints or a multicanonical 
distribution function with constraints. 

Method as claimed in claim 1 or 2, wherein the numerical sampling algorithm of at least 
one iteration generates correlated states xjj , wherein the sampling distribution function 
pj(x) preferably has a maximum that biases sampling into a sub-region of the state 
space, wherein the sampling distribution function pj(x) is preferably fitted such that the 
maximum is in a region where the product pijO{xu) has a maximum, and wherein the 
sampling preferably starts from a state close to the maximum of the sampling distribution 
function pj(x). 

Method as claimed in one of claims 4 to 6, wherein the fitting of pj(x) is done by selecting 
states x u for which the product puOfaj) has extreme values and by using the selected 
states xu to define the region in which pj{x) has extreme values. 

Method as claimed in one of claims 1 to 6, wherein parameters of the sampling distribution 
function pj{x) of at least one iteration are determined by a linear least square fit of the 
logarithm of the un-normalized sampling distribution function pj{x) to the logarithm of the 
product jfcOfat) 

Method as claimed in one of claims 1 to 8, wherein the normalization constant of the 
sampling distribution function pj(x) of at least one iteration is estimated from the sampled 
states x u and their weighting factors pij. 

Method as claimed in one of claims 1 to 9, wherein at least three iterations are done. 

Method as claimed in one of claims 1 to 10, wherein the function 0(x) is a function of a 
set of at least two functions ®^a(x)Ch(x\ . *On„(x)} 9 preferably including at least one 
of the following functions 
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at least one property for which at least one estimate is derived in the analysis of the fifth 
step, 

at least one function that is large for states that must be sampled to ensure transitions 
between important regions, 

the inverse of the probability distribution function of at least one property of the system, 
and 

the inverse of the probability distribution of the negative logarithm of the distribution 
function of the system . 

A computer software product for sampling a state space by iteratively generating states 
xis and their weighting factors pb , wherein the index / is the iteration parameter and the 
index f distinguishes different states xu generated by an iteration /, said product is 
characterized by a computer-readable medium in which program instructions are stored, 
which instructions, when read by a computer, enable the computer to 
select an initial sampling distribution function /»(x), 
to execute a fifth step for performing an analysis 
and to execute an iteration procedure including 

a second step for generating Nj states xj,i by a numerical sampling algorithm and 
a fourth step for testing at least one criterion to decide whether to continue the iteration 
procedure or to stop the iteration procedure and to go to a fifth step in order to perform the 
analysis using the simulated data, 
characterized in that the iteration further includes 

a third step determining weighting factors fat for states xij generated so far by using 
sampling distribution functions p { (x) determined so far and 

a fitting step for determining a sampling distribution function pj{x) for the next iteration by 
fitting pj(x) to ptjOfaj) for states x>,t generated so far, wherein 0(?a.t) is a function, 
respectively a property, of the states xm . 
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(57) Abstract: The present invention is an incremental 
umbrella sampling method to improve the performance 
of established sampling methods. It is sampling the 
state space by iteratively generating states xiJ and their 
weighting factors represented by Formula (a) by fitting the 
sampling distribution function p/x) of the next iteration 
to at least one weighted property of the already sampled 
states. This means that p/x) is fitted to the product 
represented by Formula (b), in which the Formula (a) 
isthe weighting factor and 0( xU ) is a function respectively 
a property of the states x u The number of states x i4 
and the number of weighting factors (see Formula (a)) 
is incremented with each iteration. In order to have a 
consistent set of weighting factors (see Formula (a)), the 
weighting factors are recalculated in each iteration for all, 
respectively for a set of selected, states. By fitting p/x) in 
the state space it is possible to use all the information of 
Formula (a) andOfc,,) for the states x ht generated so far. 
The fitting step allows to use different fitting strategies. 
For example the fitting can bias the sampling away from 
areas where intensive sampling has been done in the 
preceding iterations, or the sampling can be directed 
along local gradients respectively towards local minima or 
maxima of one or several weighted properties. In each of 
the iterations, the sampling distribution function is fitted in 
a way to improve the overall sampling of the state space. 
The method supports multi -objective optimisations. State 
space integrals can be solved. It reduces the probability that the system is trapped. The invention is general. It can be used with 
different sampling methods, in particular with Monte Carlo sampling, Metropolis Monte Carlo sampling, or dynamic simulations. 
It can be combined with the concepts of simulated annealing and multicanonical sampling. It provides a general framework that 
can be adapted to the system and the observables of interest. 
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